home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / skyline / example / dsetup.f next >
Encoding:
Text File  |  1994-08-02  |  6.4 KB  |  221 lines
  1.     subroutine dsetup( n, V, pd, nr, normV, full, band, b, job )
  2.     integer           n, job, band, pd(*), nr(*)
  3.     double precision  V(*), b(*)
  4.     double precision  normV
  5.     logical           full
  6.  
  7. c       ... Local Variables ...
  8.     integer          iv
  9.     double precision nr_sum, nr_mean, nr_std, nr_var
  10.     double precision normcol
  11.  
  12. *------------------------------------------------------------------------------
  13. *    DSETUP sets up the array V corresponding to the matrix A. The matrix
  14. *              is stored in a packed fishbone fashion as described below.
  15. *              It also sets up the rigt hand side.
  16. *
  17. *------------------------------------------------------------------------------
  18. *
  19. *    Input:
  20. *
  21. *    n :   INTEGER
  22. *             size of the problem
  23. *
  24. *       full: LOGICAL
  25. *             if .TRUE.  =>  generate a full n-by-n matrix.
  26. *             if .FALSE. =>  generate a sparse matrix.
  27. *
  28. *       band: INTEGER
  29. *             Indicates what sort of skyline to generate and the size
  30. *             of the band compared to the dimension of the problem.
  31. *
  32. *             band = 10  =>  generate a uniform banded matrix, where the
  33. *                            nonzero elements are 10% of the dimension of
  34. *                            the problem.
  35. *
  36. *             band <= 0   
  37. *                 or     =>  generate random band.
  38. *             band > 100
  39. *
  40. *             If full = .TRUE., band is ignored.
  41. *
  42. *    job : INTEGER
  43. *             Indicates the data to generate.
  44. *             job = 1  =>  generate matrix V only.
  45. *             job = 2  =>  generate rhs b only.
  46. *             job = 3  =>  generate both matrix V and rhs b.
  47. *
  48. *
  49. *       Output:
  50. *
  51. *       b : the right hand side
  52. *
  53. *    V : array where the elements are stored in a packed fashion following
  54. *        a fishbone design. That way lines below the diagonal, and columns 
  55. *        above the diagonal are stored with stride 1.
  56. *        The matrix ma has a symmetrical kind of banded-like shape. But
  57. *        the values of the matrix are not symmetrical.
  58. *
  59. *       pd : an integer array containing pointers to the diagonal elements
  60. *            of the matrix A in the array V.
  61. *       
  62. *       nr : an integer array containing the number of nonzero elements in
  63. *            a row or column up to the diagonal.
  64. *       
  65. *
  66. * V                                A
  67. * ------------------------------   ------------------------------------
  68. * |  1 | 3 | 7 | * | * | * | * |   | 11 | 12 | 13 |  * |  * |  * |  * |
  69. * |----|   |   |   |   |   |   |   |----|    |    |    |    |    |    |
  70. * |  2 | 4 | 8 |12 |18 | * | * |   | 21 | 22 | 23 | 24 | 25 |  * |  * |
  71. * |--------|   |   |   |   |   |   |---------|    |    |    |    |    |
  72. * |  5   6 | 9 |13 |19 | * |31 |   | 31 | 32 | 33 | 34 | 35 |  * | 37 |
  73. * |------------|   |   |   |   |   |--------------|    |    |    |    |
  74. * |  *  10  11 |14 |20 |24 |32 |   |  *   42   43 | 44 | 45 | 46 | 47 |
  75. * |----------------|   |   |   |   |-------------------|    |    |    |
  76. * |  *  15  16  17 |21 |25 |33 |   |  *   52   53   54 | 55 | 56 | 57 |
  77. * |--------------------|   |   |   |------------------------|    |    |
  78. * |  *   *   *  22  23 |26 |34 |   |  *    *    *   64   65 | 66 | 67 |
  79. * |------------------------|   |   |-----------------------------|    |
  80. * |  *   *  27  28  29  30 |35 |   |  *    *   73   74   75   76 | 77 |
  81. * ------------------------------   ------------------------------------
  82. *
  83. *
  84. *    pd : array of diagonal indexes
  85. *        (e.g. 1, 4, 9, 14, 21, 26, 35)
  86. *
  87. *
  88. *       For a full matrix:
  89. *
  90. *        ------------------------------------
  91. *         |  1 |  3 |  7 | 13 | 21 | 31 | 43 | 
  92. *         | 11 | 12 | 13 | 14 | 15 | 16 | 17 |
  93. *         |---------|    |    |    |    |    |
  94. *         |  2 |  4 |  8 | 14 | 22 | 32 | 44 |
  95. *         | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
  96. *         |---------------    |    |    |    |
  97. *         |  5    6 |  9 | 15 | 23 | 33 | 45 |
  98. *         | 31   32 | 33 | 34 | 35 | 36 | 37 |
  99. *         |--------------------    |    |    |
  100. *         | 10   11   12 | 16 | 24 | 34 | 46 |
  101. *         | 41   42   43 | 44 | 45 | 46 | 47 |
  102. *         |-------------------------    |    |
  103. *         | 17   18   19   20 | 25 | 35 | 47 |
  104. *         | 51   52   53   54 | 55 | 56 | 57 |
  105. *         |------------------------------    |
  106. *         | 26   27   28   29   30 | 36 | 48 |
  107. *         | 61   62   63   64   65 | 66 | 67 |
  108. *         |----------------------------------|
  109. *         | 37   38   39   40   41   42 | 49 |
  110. *         | 71   72   73   74   75   76 | 77 |
  111. *         ------------------------------------
  112. *
  113. *    pd = ( 1, 4, 9, 16, 25, 36, 49 )
  114. *
  115. *------------------------------------------------------------------------------
  116.  
  117.  
  118. c   ... Generate PD ...
  119.  
  120.     init = 1325
  121.  
  122. c   ... Full Matrix ...
  123.     if ( full ) then
  124.        do i = 1, n
  125.           pd(i) = i*i
  126.        end do
  127.  
  128.         else
  129.        pd(1) = 1
  130.  
  131.            if ( band .LE. 0  .OR.  band .GT. 100 ) then
  132.  
  133. c         ... Random bandwidth ...
  134.           do i = 2, n
  135.              pd(i) = pd(i-1) + 2*mod( 100*ran(init), i ) + 1
  136.               end do
  137.  
  138.            else
  139.  
  140. c         ... Fixed bandwidth ...
  141.               nband = n*band / 100
  142.               do i = 2, n
  143.                  pd(i) = pd(i-1) + min( nband + 1, 2*(i-1) + 1 )
  144.               end do
  145.            end if
  146.  
  147.     end if
  148.  
  149.  
  150. c   ... Calculate NR (number of nonzero elements in a row up to      ...
  151. c   ... the diagonal), its average value and the standard deviation. ...
  152.  
  153.     nr_sum = 0
  154.     nr(1)  = 0
  155.     do i = 2, n
  156.        nr(i)  = ( pd(i) - pd(i-1) ) / 2
  157.        nr_sum = nr_sum + nr(i)
  158.     end do
  159.     nr_mean = nr_sum / float(n-1)
  160.  
  161.     nr_sum = 0
  162.     do i = 1, n
  163.        nr_sum = nr_sum + ( nr(i) - nr_mean )**2
  164.         end do
  165.     nr_var = nr_sum / float(n-1)
  166.     nr_std = sqrt( nr_var )
  167.  
  168.  
  169. c   ... Generate  a matrix that will not require pivoting ...
  170.  
  171.     init = mod( 3125*init, 65536 )
  172.     V( 1 ) = ( init - 32768 ) / 16384
  173.     do j = 2, n
  174.        do i = pd(j-1)+1, pd(j)
  175.           init = mod( 3125*init, 65536 )
  176.           V( i ) = ( init - 32768 ) / 16384
  177.        end do
  178.     end do
  179.     do j = 1, n
  180.        normcol = 0.d0
  181.        do i = pd(j) - nr(j), pd(j)
  182.           normcol = normcol + abs( V( i ) )
  183.        end do
  184.        do i = j+1, n
  185.           if( nr(i) .GE. i-j ) 
  186.      &            normcol = normcol + abs( V( pd(i) - nr(i) -(i-j) ) )
  187.        end do
  188.        V( pd( j ) ) = V( pd( j ) ) + normcol + 1.d0
  189.     end do
  190.  
  191.  
  192. c  ...  Calculate the norm of A. ...
  193.  
  194.     normV = max( V( 1 ), 0.d0 )
  195.     do j = 2, n
  196.        do i = pd( j-1 ) + 1, pd( j )
  197.           normV  = max( V( i ), normV )
  198.        end do
  199.     end do
  200.  
  201.  
  202. c   ... Generate the right hand side ...
  203.  
  204.     do i = 1, n
  205.        b(i) = 0.d0
  206.     end do
  207.  
  208.     do j = 1, n
  209.        do i = pd(j) - nr(j), pd(j)
  210.           b( j - pd(j) + i ) = b( j - pd(j) + i ) + V(i)
  211.        end do
  212.        do i = j+1, n
  213.           if( nr(i) .GE. i-j )
  214.      &              b(i) = b(i) + V( pd(i) - nr(i) -(i-j) )
  215.        end do
  216.     end do
  217.  
  218.  
  219.     return
  220.     end
  221.